!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine model_variables (ro, mx, my, mz, bx, by, bz, en, gx, gy, po, x, y, nx, ny)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny

   real(8), dimension(nx), intent(in) :: x
   real(8), dimension(ny), intent(in) :: y
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy, po

   real(8), dimension(nx, ny) :: te, pr
   real(8), dimension(ny_total_finest) :: te_total_y

   integer(4) :: i, j, j_total

   real(8) :: x0, y0, r2, v_amp, c_x
   real(8) :: par_b, par_i, r0, mag_e, beta, g0, te_chr, te_cor, h_cor, w_tr, zeroy, dy

!----------------------------------------------------------------------------------------------------------------------------------|
!  density, temperature,
!----------------------------------------------------------------------------------------------------------------------------------|
   g0 = -1.0d0
   te_chr = 1.0d0
   te_cor = 100.0d0
   h_cor = 10.0d0
   w_tr = 1.0d0

   if (.not. allocated (y_total)) allocate (y_total(ny_total_finest))
   if (.not. allocated (ro_total_y)) allocate (ro_total_y(ny_total_finest))

   if (ro_total_y_flag) then

      zeroy = (nint ((ny_total - 2 * ng) * zero_posy / rangey) + ng) * 2 ** (n_levels - 1) + 1
      dy = rangey / dble(ny_total - 2 * ng) / 2.0d0 ** (n_levels - 1)

      y_total(zeroy) = dy * 0.5d0

      do j = zeroy + 1, ny_total_finest
         y_total(j) = y_total(j - 1) + dy
      enddo
      do j = ny_total_finest - 1, 1, -1
         y_total(j) = y_total(j + 1) - dy
      enddo

      do j = 1, ny_total_finest
         te_total_y(j) = te_chr + (te_cor - te_chr) / 2.0d0 * (tanh ((y_total(j) - h_cor) / w_tr) + 1.0d0 + eps)

         if (j .eq. 1) then
            ro_total_y(j) = 1.0d3
         else
            ro_total_y(j) = (2.0d0 * te_total_y(j - 1) + (y_total(j) - y_total(j - 1)) * g0) /                                     &
               (2.0d0 * te_total_y(j) - (y_total(j) - y_total(j - 1)) * g0) * ro_total_y(j - 1)
         endif
      enddo

      ro_total_y_flag = .false.

   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  magnetic flux function
!----------------------------------------------------------------------------------------------------------------------------------|
   beta = 0.5d0
   r0 = 0.5d0
   par_b = 1.2d2
   par_i = par_b * 1.2d0
   mag_e = 0.5d0
   v_amp = 1.0d0
   c_x = 40.0d0

   if (.not. allocated (x_total)) then
      allocate (x_total(1))
      x_total = c_x
   endif

   do j = 1, ny
   do i = 1, nx

      gx(i, j) = 0.0d0
      gy(i, j) = g0
      bx(i, j) = 0.0d0
      by(i, j) = 0.0d0
      bz(i, j) = 0.0d0

!----------------------------------------------------------------------------------------------------------------------------------|
!     background field
!----------------------------------------------------------------------------------------------------------------------------------|
      x0 = -x_total(1) 
      y0 = -10.0d0
      r2 = (x(i) - x0) ** 2 + (y(j) - y0) ** 2
      bx(i, j) = bx(i, j) + par_b * (y(j) - y0) / r2
      by(i, j) = by(i, j) - par_b * (x(i) - x0) / r2

      x0 = x_total(1)
      y0 = -10.0d0
      r2 = (x(i) - x0) ** 2 + (y(j) - y0) ** 2
      bx(i, j) = bx(i, j) + par_b * (y(j) - y0) / r2
      by(i, j) = by(i, j) - par_b * (x(i) - x0) / r2

!----------------------------------------------------------------------------------------------------------------------------------|
!     image current
!----------------------------------------------------------------------------------------------------------------------------------|
      x0 = 0.0d0
      y0 = -20.0d0
      r2 = (x(i) - x0) ** 2 + (y(j) - y0) ** 2
      bx(i, j) = bx(i, j) - par_i * (y(j) - y0) / r2
      by(i, j) = by(i, j) + par_i * (x(i) - x0) / r2

!----------------------------------------------------------------------------------------------------------------------------------|
!     other variables
!----------------------------------------------------------------------------------------------------------------------------------|
      te(i, j) = te_chr + (te_cor - te_chr) / 2.0d0 * (tanh ((y(j) - h_cor) / w_tr) + 1.0d0 + eps)

      do j_total = 1, ny_total_finest
         if (y_total(j_total) .ge. y(j)) exit
      enddo

      if (j_total .eq. 1) then
         ro(i, j) = ro_total_y(j_total)
      else
         ro(i, j) = (ro_total_y(j_total - 1) * (y_total(j_total) - y(j)) + ro_total_y(j_total) * (y(j) - y_total(j_total - 1))) /  &
            (y_total(j_total) - y_total(j_total - 1))
      endif

      pr(i, j) = ro(i, j) * te(i, j)

      mx(i, j) = 0.0d0
      my(i, j) = 0.0d0
      if (y(j) .le. 0.0d0 .and. abs (x(i)) .le. 2.0d0 * c_x) mx(i, j) = -v_amp * sin (x(i) * pi / c_x / 2.0d0) * ro(i, j)
!      if (abs (x(i)) .le. 2.0d0 * c_x .and. abs (y(j)) .le. 20.0d0) mx(i, j) = -v_amp * sin (x(i) * pi / c_x / 2.0d0) * cos (y(j) * pi / 40.0d0) * ro(i, j)
      mz(i, j) = 0.0d0

      en(i, j) = (mx(i, j) * mx(i, j) + my(i, j) * my(i, j) + mz(i, j) * mz(i, j)) / ro(i, j) / 2.0d0 +                            &
         (bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)) / 2.0d0 + pr(i, j) / gmm1

      po(i, j) = 0.0d0

   enddo
   enddo
   x_total = c_x - v_amp * sin (x_total * pi / c_x / 2.0d0) * 1.0d-4

!----------------------------------------------------------------------------------------------------------------------------------|
end subroutine model_variables
